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Abstract. For decades, the genus Pelomedusa has been thought of as containing only a single species that is widely distrib¬ 
uted in sub-Saharan Africa, the Arabian Peninsula, and Madagascar. Elowever, ten species and up to five candidate spe¬ 
cies have been recognized recently. For these taxa, ecological information is scarce and distribution ranges are often only 
incompletely known. In this study, ecological niche spaces were estimated for eight Pelomedusa species, using a maximum 
entropy algorithm and geographical, climatic, and topographic information. Using the obtained models and Schoener s D 
and Elellinger’s I indices, areas of sympatry between species were estimated and candidate areas for possibly overlooked 
taxa were identified. Furthermore, differences in niche spaces between species were assessed. The ecological niche models 
revealed geographical spaces that match the predicted niche of the eight examined species. However, for P. subrufa sensu 
stricto, the distribution range of the introduced populations in Madagascar was not predicted. Our models, together with 
statistical analyses, indicate that the studied Pelomedusa species occupy different geographical spaces. Areas of sympatry 
with significant niche overlap were identified between some species pairs. Yet, Pelomedusa species are generally allopatri- 
cally or parapatrically distributed. A minimum of five candidate regions possibly harbouring unknown taxa were identi¬ 
fied. Our results provide a better understanding of the distribution ranges of Pelomedusa species and, thus, can serve as a 
basis for conservation planning. 
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Introduction 

Complexes of cryptic species require special consideration 
in conservation planning because unrecognised species can 
be rare, endangered, and may therefore need special con¬ 
servation strategies (Bickford et al. 2007). Assessing dis¬ 
tribution patterns of cryptic taxa is often difficult because 
many locality records can be misidentified. Helmeted ter¬ 
rapins (genus Pelomedusa) represent an excellent example 
for this: For a long time, these terrapins were identified as 
a single widespread species, Pelomedusa subrufa (Bonna- 
terre, 1789) sensu lato. It was thought to be distributed in 
a vast range, encompassing sub-Saharan Africa, the south¬ 
western Arabian Peninsula, and Madagascar (Iverson 
1992, Gasperetti et al. 1993). However, recent studies have 
demonstrated that helmeted terrapins represent a diverse 
species complex (VARCAS-RAMfREZ et al. 2010, Wong et al. 
2010, Fritz et al. 2011, 2014, Petzold et al. 2014, Nagy et 
al. 2015) with currently 10 named species and up to another 


five unnamed candidate species in continental Africa and 
Arabia (Fritz et al. 2014, Petzold et al. 2014, Nagy et al. 
2015). The helmeted terrapins from Madagascar are now 
understood as introduced populations of P. subrufa sensu 
stricto (VARGAS-RAMfREZ et al. 2010, Wong et al. 2010, 
Petzold et al. 2014). However, without genetic verifica¬ 
tion, many of the old records of helmeted terrapins cannot 
be identified with a certain species, which is why the distri¬ 
bution ranges of many taxa are only incompletely known. 

Based on known locality data, distribution modelling 
can improve this situation. Among other biogeographi- 
cal analyses (Raxworthy et al. 2007, Rissler & Apodaca 
2007, Nakazato et al. 2010), distribution modelling relat¬ 
ing verified locality and environmental data has generally 
been applied to identify regions where rare or undiscov¬ 
ered species may occur (Anderson et al. 2002, Ferrier 
et al. 2002, Raxworthy et al. 2003, 2007, Siqueira et al. 
2009), to recognize environmental differences among spe¬ 
cies (Kozak et al. 2008), to identify regions of high ende- 
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micity and conservation value (Rissler et al. 2006), and 
to support conservation planning (e.g., Araujo & Wil¬ 
liams 2000, Perrier et al. 2002). Based on the premise 
that the observed distribution of a species provides useful 
information as to its environmental requirements (Pear¬ 
son & Dawson 2003), this approach associates known dis¬ 
tributional information with a set of environmental vari¬ 
ables to characterize the range of habitat conditions (the 
realized niche sensu Hutchinson 1957). Thus, ecological 
niche modelling (ENM) utilizes environmental conditions 
at known presence localities of a target species to assem¬ 
ble a model of its ecological tolerances that can be used to 
predict the geographical space that fulfils the same require¬ 
ments (Ackerly 2003, Pearson et al. 2007). However, the 
results from ENM must be interpreted cautiously because 
factors like geographical barriers, biotic interactions, and 
demographic history are often not considered, implying 
that species rarely occupy all areas with suitable environ¬ 
ments (Anderson et al. 2002, Pearson & Dawson 2003, 
Araujo & Pearson 2005, Phillips et al. 2006). Regardless 
of these limitations, distribution models have been shown 
to provide highly informative biogeographical information 
(e.g., Pleishman et al. 2002, Bourg et al. 2005, Raxwor- 
THY et al. 2007), especially if researchers adopt a conserva¬ 
tive interpretation of model results as areas similar to those 
from which a species is known (Pearson et al. 2007). 

Based on genetically verified locality data for Pelomedu¬ 
sa species (Petzold et al. 2014), we perform in the present 
paper ENM and climatic niche analyses (i) to predict the 
ecological niches and geographical spaces for several spe¬ 
cies of helmeted terrapins, (ii) to identify areas of sympa- 
try, and (in) to assess ecological variation among species. 
However, we excluded from our modelling approach the 
unnamed candidate species and two named species (P ge- 
hafie, P. somalica) for which only one or two collection 
sites are known (Petzold et al. 2014, Nagy et al. 2015) due 
to the sensitivity of ecological niche model approaches to 
small sample sizes (Pearson et al. 2007, Wisz et al. 2008). 
As an additional aim, we try to identify areas with further, 
as yet possibly overlooked Pelomedusa species. 

Materials and methods 

Estimation of climatic niche 

The climatic niche was analysed in two ways. In an ENM 
approach, the climate niche was projected across a geo¬ 
graphical region using the georeferenced locality data for 
each Pelomedusa species from Petzold et al. (2014). How¬ 
ever, locality data for the introduced Malagasy popula¬ 
tions of P. subrufa sensu stricto were excluded. In addition, 
a point-based analysis was conducted utilizing univari¬ 
ate and multivariate statistics for specific conditions using 
climate variables and altitude alone. Bioclimatic variables 
and altitude were taken from WorldClim version 1.4 with 
a resolution of 30 arcseconds (Hijmans et al. 2005). In or¬ 
der to reduce effects of overfitting and multi-collinearity 
of predictors (Heikkinen et al. 2006), correlated variables 


were identified using R v.2.9.2. Por doing so, climatic and 
altitude data across Africa and the Arabian Peninsula were 
extracted for 10^ randomly generated points using the ‘Cre¬ 
ate Random Points command in ArcGIS v.10.2. (ESRI Spa¬ 
tial Analyst), and a Pearson correlation coefficient was es¬ 
timated for each pair of variables. Due to the expanse of 
the area of study, this procedure was repeated three times. 
For the pair of variables that exceeded a threshold of r = 

O. 75 (Kalkvik et al. 2011), the variable that was biologi¬ 
cally most meaningful was integrated. Based on a the re¬ 
sults of the three repetitions, the following eleven variables 
were included in the analyses: (1) annual mean tempera¬ 
ture, (2) mean diurnal range, (3) isothermality, (4) max¬ 
imum temperature of the warmest month, (5) minimum 
temperature of the coldest month, (6) annual precipitation, 
(7) precipitation of the wettest quarter, (8) precipitation of 
the driest quarter, (9) precipitation of the warmest quarter, 
(10) precipitation of the coldest quarter, and (11) altitude. 

Niche modelling, overlap, and areas of distribution 

Niche modelling was used to identify areas with climatic 
and topographic conditions resembling those regions for 
which verified records of the different Pelomedusa species 
exist. Models were created using MAXENT v.3.3.3k, which 
applies a maximum entropy algorithm to estimate dis¬ 
tributions based on niche characteristics (Phillips et al. 
2006, Phillips & DudIk 2008). Maximum entropy is use¬ 
ful when absence data is lacking and performs better than 
other distribution modelling algorithms when the number 
of data points is limited (Elith et al. 2006, Pearson et 
al. 2007, Wisz et al. 2008). Our niche modelling includ¬ 
ed Africa, Madagascar, and the Arabian Peninsula, i.e., the 
whole distribution range of all Pelomedusa species. Only 
a few genetically verified occurrence points are available 
for Pelomedusa barbata, P. kobe, P. neumanni, P. olivacea, 

P. schweinfurthi, and P. variabilis (Petzold et al. 2014). To 
increase sample size for these species, a minimum con¬ 
vex polygon (MCP) was created using the genetically ver¬ 
ified records from Petzold et al. (2014). This sampling 
was enriched with additional Pelomedusa records within 
each polygon, which were allocated to the respective spe¬ 
cies (Supplementary figs S1-1-S2). In order to increase the 
number of sites for P. barbata and P. schweinfurthi, records 
immediately adjacent to the polygons were also includ¬ 
ed (Table 1). Additional locality data were obtained from 
HerpNET (http://www.herpnet.org) and EMYSystem 
(http://emys.geo.orst.edu/; Supplementary table Si). To as¬ 
sess the accuracy of each model and due to the low number 
of localities, the jackknife data partitioning approach 
(Pearson et al. 2007) was used. As described by these au¬ 
thors, independent MAXENT models were run with one 
of the localities excluded in each case. For each model, the 
lowest presence threshold was used and the ability to pre¬ 
dict the excluded locality assessed. For doing so, the pro¬ 
gram PVALUECOMPUTE.EXE (Pearson et al. 2007) 
was used, and a p-value was calculated for each species. 
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Table 1. Model contributions of the bioclimatic variables and altitude for the Pelomedusa species (percentages) and results of model 
evaluation. Bioclimatic variables: ALT - altitude; AMT - annual mean temperature; AP - annual precipitation; IT - isothermality; 
MDR - mean diurnal range; MTCM - minimum temperature of the coldest month; MTWM - maximum temperature of the warmest 
month; PCQ - precipitation of the coldest quarter; PDQ - precipitation of the driest quarter; PWaQ - precipitation of the warmest 
quarter; PWQ - precipitation of the wettest quarter, n - number of records; O.r. - omission rate; AUC - area under the curve of 
ROC plots. 


Rank 

barbata 


galeata 


kobe 


neumanni 

olivacea 


schweinfurthi 

subrufa 


variabilis 

1 

MTCM 

26.2 

PDQ 

34.7 

IT 

22.6 

IT 

57.1 

PCQ 

21.0 

IT 

31.1 

MTCM 

34.3 

MTCM 

73.5 

2 

MDR 

20.1 

IT 

18.4 

ALT 

22.3 

ALT 

22.0 

PWQ 

19.6 

PCQ 

18.8 

ALT 

27.1 

MDR 

15.8 

3 

AP 

19.3 

AMT 

16.9 

PDQ 

16.8 

MTCM 

5.9 

PDQ 

14.6 

PWaQ 

16.1 

MTWM 11.8 

ALT 

4.5 

4 

PWQ 

17.0 

MTCM 

11.1 

PCQ 

14.2 

AP 

2.7 

AMT 

11.8 

MTCM 

14.0 

PCQ 

6.0 

AP 

3.9 

5 

ALT 

13.1 

ALT 

7.0 

MDR 

10.1 

MDR 

2.5 

MDR 

11.6 

PWQ 

7.6 

IT 

5.5 

AMT 

0.5 

6 

PWaQ 

1.6 

PWaQ 

4.1 

MTWM 

8.0 

AMT 

2.5 

PWaQ 

10.0 

MDR 

4.9 

PWaQ 

4.3 

IT 

0.4 

7 

PDQ 

0.9 

MDR 

3.4 

AP 

3.4 

MTWM 

2.3 

ALT 

5.1 

AP 

4.9 

MDR 

3.5 

PDQ 

0.4 

8 

AMT 

0.7 

PWQ 

1.6 

PWaQ 

1.3 

PWaQ 

1.8 

AP 

2.6 

ALT 

1.3 

PDQ 

2.5 

PWQ 

0.4 

9 

PCQ 

0.7 

MTWM 

1.0 

MTCM 

1.1 

PDQ 

1.7 

IT 

1.5 

PDQ 

1.1 

AMT 

2.3 

PCQ 

0.3 

10 

IT 

0.3 

AP 

0.9 

PWQ 

0.1 

PCQ 

1.2 

MTCM 

1.2 

MTWM 

0.1 

PWQ 

1.7 

MTWM 

0.2 

11 

MTWM 

0.1 

PCQ 

0.8 

AMT 

0.1 

PWQ 

0.3 

MTWM 

1.0 

AMT 

0.0 

AP 

1.0 

PWaQ 

0.1 

n 

11 


23 


10 


13 


18 


10 


16 


12 



1 of 11 


2 of 23 


1 of 10 


2 of 13 


2 of 18 


1 of 10 


1 of 16 


1 of 12 



p<0.001 


p<0.001 


p<0.001 


p<0.001 


p<0.001 


p<0.001 


p<0.001 


p<0.001 


AUC 

0.991 


0.997 


0.984 


0.990 


0.990 


0.992 


0.955 


0.991 



Additionally, because all included species were represented 
by more than ten locality records, a model evaluation was 
performed using the area under the curve (AUC) of the 
receiver-operating characteristic (ROC) plot. For this, the 
occurrence data were randomly partitioned into training 
(75%) and test (25%) data sets. The results of each model 
were summarized as averages of 100 bootstrap replicates. 
AU analyses were run with a convergence threshold of 1.0 x 
10^ with 1,000 iterations. The regularization multiplier was 
set at 0.5 to obtain a more localized output distribution that 
will be a closer fit to the given presence records (Rivera 
et al. 2011). For each species, the minimum training pres¬ 
ence was used as the threshold for generating binary pre¬ 
diction maps of suitable climatic niche. A jackknife resam¬ 
pling test for evaluating the importance of each variable 
was performed, and to examine how each variable affected 
the MAXENT prediction, individual response curves were 
examined. Disjunct areas of major overprediction were ex¬ 
cluded from further analyses in this context (but see be¬ 
low). Niche overlaps of lineages were estimated using two 
different measures of similarity. Schooner’s D and the I sta¬ 
tistic, based on Ffellinger distance. Each metric is calcu¬ 
lated by comparing the estimates of habitat suitability from 
the ENMs generated by MAXENT for each grid cell of the 
study area after normalizing each ENM, so that all suit¬ 
ability scores sum up to 1. Both similarity tests range from 
o (no overlap between predicted environmental toleranc¬ 
es of species) to 1 (all grid cells are estimated to be equal¬ 
ly suitable for both species). The ecological interpretation 
of D implies that the suitability scores are proportional to 
species abundance, while I simply treats the two ENMs as 
probability distributions (Warren et al. 2010). Further¬ 


more, range overlap between species was calculated using 
the average of all minimum training presences as a thresh¬ 
old. This metric and niche overlap were computed using 
ENMTOOLS v.1.4.3 (Warren et al. 2008, 2010; available at 
http://enmtools.blogspot.com). In addition, the geographi¬ 
cal distribution of each species was estimated by means of 
the extent of occurrence (EOO) and the area of occupan¬ 
cy (AOO). These estimators are used in some criteria for 
determining the threat category of a species (lUCN 2012). 
The EOO for each species was calculated in km^ from the 
MCPs using the GeoCAT portal (Bachman et al. 2011; 
http://geocat.kew.org/) and, in addition, using ArcGIS 
from the geographical space that fulfils the predicted eco¬ 
logical niche. The AOO was calculated using the GeoCAT 
portal. The EOO is defined as “the area contained within 
the shortest continuous imaginary boundary, which can be 
drawn to encompass aU the known, inferred or projected 
sites of present occurrence of a taxon”, and the AOO “as the 
area within its extent of occurrence’, which is occupied by a 
taxon” (lUCN 2012), for both, excluding cases of vagrancy. 

Identification of candidate areas for possibly unknown 
additional Pelomedusa species 

For some species, ecological niche modelling has guided 
the discovery of isolated areas of environmental suitability 
that are not actually occupied by the species being mod¬ 
elled (Raxworthy et al. 2007). Species rarely occupy all 
regions with suitable environments due to historical and 
ecological constraints like geographical barriers or biot¬ 
ic interactions (e.g., interspecific competition, predation 
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pressure or dispersal abilities). Such factors, which can 
play a key role in shaping the geographical ranges of spe¬ 
cies (Peterson 2001, Svenning & Skov 2004, Araujo & 
Pearson 2005, Raxworthy et al. 2007), are not consid¬ 
ered in our modelling approach. However, disjunct over¬ 
predicted regions can be valuable indicators for unrecog¬ 
nised distribution ranges or even overlooked taxa. Rax¬ 
worthy et al. (2003,2007) have pointed out that especially 
in poorly studied regions with many local endemics, iso¬ 
lated areas of overprediction can correspond to the ranges 
of unknown taxa. Here we adopt this approach as a means 
to identify candidate areas where possibly unknown spe¬ 
cies of Pelomedusa might occur. For doing so, we com¬ 
pare our disjunct areas of overprediction with published 
locality records for helmeted terrapins from the HerpNET 
data portal (http://www.herpnet.org) and the world turtle 
database of the EMYSystem data portal (http://emys.geo. 
orst.edu). All used records are shown in Supplementary 
figure S2. We hypothesize that isolated areas of overpre¬ 
diction potentially harbour additional candidate species 
of Pelomedusa, especially if records of helmeted terrapins 
exist from these areas. Moreover, records of helmeted ter¬ 
rapins both outside the contiguous and isolated regions 
predicted for the studied Pelomedusa species could corre¬ 
spond to additional unrecognised taxa. 

Point-based analysis 

To establish whether significant differences exist among 
the climatic spaces of different species, univariate (ANO- 
VA) and multivariate analyses of variance (MANOVA) and 
parametric discriminant analysis (PDA) were performed 
using occurrence points. For comparing the multidimen¬ 
sional climatic values between records, altitude and climat¬ 
ic values were extracted for each site using ArcGIS. A one¬ 
way analysis of variance (ANOVA) was used on each of the 
10 environmental variables and altitude to test for signifi¬ 
cant differences between species. Subsequently, a MANO¬ 
VA was calculated using all variables simultaneously to as¬ 
sess whether significant overall differences exist among the 
ecological spaces occupied by the individual species. Fur¬ 
thermore, the relative contributions of environmental vari¬ 
ables and altitude to the discrimination of the species were 
quantified by a PDA. Divergences were visualized using 
the centroid for each lineage. All sampling locations repre¬ 
senting explicit climatic and altitude data were plotted on 
the main discriminant dimensions. All statistical analyses 
were performed in SPSS v.19.0. 

Results 

Estimation of climatic niche and niche modelling, 
overlap, and areas of distribution 

The jackknife data partitioning validation revealed high 
success rates (i.e., low omission rates) and showed that all 
models were significantly better than random expectations 


(Table 1). Furthermore, the AUC values for the test data 
set were robust and revealed high specificity for each mod¬ 
el. The minimum probability of occurrence for training 
points ranged from 0.017 for Pelomedusa galeata to 0.381 
for P. subrufa (Supplementary figs S3-S10). In the obtained 
models, the predicted ecological niche of each species was 
affected, with different strength, by virtually all environ¬ 
mental variables (Table 1). The minimum temperature of 
the coldest month made the largest contribution to P. bar- 
bata, P. subrufa, and P. variabilis, isothermality was most 
important for P. kobe, P. neumanni, and P. schweinfurthi, 
precipitation of the coldest quarter was most important 
for P. olivacea, and precipitation of the driest quarter for 
P. galeata (Table 1). 

For characterizing the geographical ranges fulfilling the 
predicted ecological niches of the studied Pelomedusa spe¬ 
cies (Fig. 1), the classification of sub-Saharan regions by 
Linder et al. (2012) will be followed here. The predicted 
distribution of P. olivacea ranges along the western Suda- 
nian Region. The predicted main area included a savan¬ 
nah belt from Senegal to southern Chad (Fig. lA). A ma¬ 
jor, significantly overpredicted disjunct area was revealed 
between southern Sudan and northern South Sudan (Sup¬ 
plementary fig. S3), corresponding to the eastern Sudanian 
Region. Additional but smaller disjunct overpredicted ar¬ 
eas were also indicated along the coast of Angola, in east¬ 
ern Zambia and northernmost Zimbabwe, in northern 
Mozambique, Malawi, Tanzania, Kenya, Somalia, Ethiopia, 
and in western Madagascar. For P. neumanni, the mod¬ 
el predicted two main areas (Fig. lA): one located within 
the Ethiopian and Somalian Regions, and a second one in 
the northern Zambezian Region of the eastern Democrat¬ 
ic Republic of the Congo, Uganda, Tanzania, and Kenya 
(Eig. lA). Disjunct small overpredicted areas were indicat¬ 
ed in northern and central Madagascar, northern Somalia, 
and western Angola (Supplementary fig. S4). The occur¬ 
rence of P. galeata was predicted all over the Southern Af¬ 
rica Region, extending into the southeastern corner of the 
Zambezian Region (Eig. lA). Disjunct small overpredicted 
areas were indicated along the coast of the western Magh¬ 
reb, in the southwestern Arabian Peninsula, Eritrea, Ethi¬ 
opia, Somalia, Kenya, Malawi, and southern Madagascar 
(Supplementary fig. S5). 

The predicted distribution of P. barbata extends along 
the southern coast of the Arabian Peninsula from Saudi 
Arabia through Yemen to Oman (Fig. iB). Further disjunct 
overpredicted areas were revealed on the southern Red 
Sea coast of Africa, the coastal regions of the Horn of Af¬ 
rica, and inland in the border region of Kenya and Ethio¬ 
pia (Supplementary fig. S6). The occurrence of P. schwein¬ 
furthi was predicted for a large region encompassing parts 
of eastern and southern Ethiopia, Sudan, South Sudan, the 
northeastern Democratic Republic of the Congo, north¬ 
ern Uganda, and northern Kenya and adjacent Somalia 
(Eig. iB). This range corresponds to the crossroads of sev¬ 
eral biogeographical regions of Linder et al. (2012). Small 
disjunct overpredicted areas for P. schweinfurthi were 
suggested for several regions in West Africa and the cen- 
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tral Democratic Republic of the Congo (Supplementary 
fig. S7). Pelomedusa subrufa was predicted for a large area 
corresponding to the Zambezian Region and the Kalahari 
Subregion, and parts of the Natal and Namib Subregions of 
the Southern African Region of Linder et al. (2012). Small 
disjunct overpredicted areas were also revealed for Chad, 
Sudan, the East African countries, the southwestern Ara¬ 
bian Peninsula, and inland Madagascar (Supplementary 
fig. S8). However, the species was not predicted for any of 
the regions in northern, western and southern Madagas¬ 
car, from where introduced populations of P. subrufa are 
known (cf the EMYSystem map and Petzold et al. 2014). 


The occurrence of P. variabilis (Pig. rC) was predicted 
for the Guinean Subregion and the coastal Congolian Sub- 
region of the Congolian Region of Linder et al. (2012). 
In addition, major suitable areas were suggested for the 
southern and northern Red Sea coasts, southern Soma¬ 
lia and adjacent Kenya and coastal Tanzania, northwest¬ 
ern Kenya, and the border region of the Congo Republics 
(Supplementary fig. S9). Por P. kobe, the northeastern part 
of the Zambezian region and parts of the Somalian Region 
were revealed as suitable (Fig. rC). Disjunct overpredicted 
areas were found in western Angola and Namibia as well as 
Yemen (Supplementary fig. Sro). 



Figure 1. Environments modelled as suitable using MAXENT for A) Pelomedusa olivacea (yellow), P. neumanni (green), and P. galeata 
(brown); B) P. barbata (black), P. schweinfurthi (purple), and P. subrufa (blue); and C) P. variabilis (orange) and P. kobe (grey). Ranges 
are based on 10 environmental variables and altitude. Minimum training presence thresholds used to separate the map into binary 
predictions (suitable: coloured, unsuitable: not coloured). Major overpredicted disjunct areas except Madagascar are not shown. Black 
(or white) dots: records used for developing the models. Scale bars = 1,000 km. Inset (left): P. galeata, photo: W. R. Branch, inset 
(right): P. subrufa, photo: A. Schleicher. 
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Table 2. Percentage of geographical overlap in niche spaces pre¬ 
dicted by MAXENT and measures of niche overlap (Schoener s 
D and Hellingers I indices) between eight species of Pelomedusa. 


Species pair 

% overlap 

I 

D 

P. barbata vs P. galeata 

0.0 

0.000 

0.000 

P. barbata vs P. kobe 

0.0 

0.000 

0.000 

P. barbata vs P. neumanni 

0.0 

0.000 

0.000 

P. barbata vs P. olivacea 

0.0 

0.000 

0.000 

P. barbata vs P. schweinfurthi 

0.0 

0.000 

0.000 

P. barbata vs P. subrufa 

0.0 

0.000 

0.000 

P. barbata vs P. variabilis 

0.0 

0.000 

0.000 

P. galeata vs P. kobe 

0.0 

0.000 

0.000 

P. galeata vs P. neumanni 

0.0 

0.000 

0.000 

P. galeata vs P. olivacea 

0.0 

0.000 

0.000 

P. galeata vs P. schweinfurthi 

0.0 

0.000 

0.000 

P. galeata vs P. subrufa 

39.6 

0.266 

0.176 

P. galeata vs P. variabilis 

0.0 

0.000 

0.000 

P. kobe vs P. neumanni 

25.3 

0.222 

0.179 

P. kobe vs P. olivacea 

0.0 

0.000 

0.000 

P. kobe vs P. schweinfurthi 

3.8 

0.031 

0.025 

P. kobe vs P. subrufa 

8.3 

0.043 

0.019 

P. kobe vs P. variabilis 

0.0 

0.000 

0.000 

P. neumanni vs P. olivacea 

0.0 

0.000 

0.000 

P. neumanni vs P. schweinfurthi 

4.4 

0.037 

0.030 

P. neumanni vs P. subrufa 

7.2 

0.042 

0.020 

P. neumanni vs P. variabilis 

0.0 

0.000 

0.000 

P. olivacea vs P. schweinfurthi 

0.0 

0.000 

0.000 

P. olivacea vs P. subrufa 

0.0 

0.000 

0.000 

P. olivacea vs P. variabilis 

1.6 

0.004 

0.002 

P. schweinfurthi vs P. subrufa 

0.0 

0.000 

0.000 

P. schweinfurthi vs P. variabilis 

0.0 

0.000 

0.000 

P. subrufa vs P. variabilis 

0.0 

0.000 

0.000 


The revealed environmental niches of the individual 
Pelomedusa species are mostly non-overlapping or at least 
little overlapping, indicating that allopatry or parapatry 
constitute the general spatial configurations (Fig. 2, Ta¬ 
ble 2). However, a large area of geographical niche over¬ 
lap was revealed in eastern South Africa for P. galeata and 
P. subrufa (39.6%; I = 0.266, D = 0.176). Moreover, signifi¬ 
cant overlap was also suggested for the East African spe¬ 
cies P. kobe and P. neumanni (25.3%; I = 0.222, D = 0.179). 
The largest extents of occurrence (EOO) were found for 
P. subrufa, P. olivacea, and P. galeata, and the smallest for 
P. barbata and P. variabilis (Table 3). Eurthermore, distinct¬ 
ly smaller areas of occupancy (AOO) were revealed for all 
species (Table 3). 

Identification of candidate areas for possibly unknown 
additional Pelomedusa species 

Seven disjunct overpredicted major areas (I-VII in Fig. 3) 
containing records for helmeted terrapins were identi- 


Table 3. Extents of occurrence (EOO) and areas of occupancy 
(AOO) in km^. MCP - minimum convex polygon. 


Species 

EOO 

MAXENT 

MCP 

AOO 

P. barbata 

320190 

66154 

44000 

P. galeata 

1435633 

858795 

92000 

P. kobe 

860929 

83956 

40000 

P. neumanni 

991562 

128206 

48000 

P. olivacea 

2283126 

1292791 

72000 

P. schweinfurthi 

906026 

180989 

40000 

P. subrufa 

2674085 

3313210 

64000 

P. variabilis 

475410 

99671 

40000 


fied. These areas often correspond well to the statistical¬ 
ly defined biogeographical regions of sub-Saharan Africa 
of Linder et al. (2012). Area I matches largely the Namib 
and southwest Angola Subregions of the Southern Afri¬ 
can Region, area II falls within the Congolian Subregion 
of the Congolian Region, and area III is located in the east¬ 
ern part of the Sudanian Region. In this area lies the only 
known collection site for candidate species B of Petzold 
et al. (2014), which represents a highly distinct genetic line¬ 
age. Candidate species B was not described as a new species 
because only one specimen is known. Area IV corresponds 
to the Ethiopian Subregion, and area V largely matches the 
Horn Subregion of the Somalian Region. Area V embraces 
the two genetically verified records of P. somalica. This spe¬ 
cies consists of two highly distinct genetic lineages record¬ 
ed from two neighbouring sites; it corresponds most likely 
to two distinct taxa (Petzold et al. 2014). Due to the pau¬ 
city of records, P. somalica was not included in the present 
modelling. Area VI lies in the eastern Somalian Region 
and the northeastern Zambezian Region, and area VII in 
the eastern inland portion of the Zambezian Region. 

Additional records for helmeted terrapins fall neither 
within any of the areas modelled for the eight Pelomedusa 
species nor in any major disjunct area of overprediction 
(Pigs 1+3). Among others, such outliers were identified for 
the Democratic Republic of the Congo. A new candidate 
species was recently described from the southeastern part 
of this country (Nagy et al. 2015). 

Point-based analysis 

One-way analysis of variance (ANOVA) revealed for all 
ten climate variables and altitude significant differences 
between species (Supplementary table S2). Furthermore, 
multivariate analyses of variance (MANOVA) indicated 
significant overall differences among the ecological spaces 
of the eight species (Wilks’ lambda = 0.0001, F = 23.083, 
p < 0.001). The two statistical tests indicate species-spe¬ 
cific ecological preferences. This is also underlined by the 
PDA results (Fig. 4). The two first canonical scores of the 
PDA explained 80.8% of environmental variation among 
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the eight studied species. For the 61.5% of the variation 
explained by the first axis, minimum temperature of the 
coldest month and isothermality were the most important 
variables. For the 19.3% of the variation explained by the 
second axis, maximum temperature of the warmest month 
and annual mean temperature were most important. The 
first axis shows a gradient from subtropical sites with 
marked cold winters to tropical sites with warmer temper¬ 
atures, and the second axis reflects a gradient from regions 
with low annual mean temperatures to regions with high 
annual mean temperatures. Another MANOVA analysis, 
using the classification function coefficients as dependent 
variables and species as categorical variables, confirmed 



differences in the environmental spaces of the eight species 
(Wilks’ lambda = 0.526, F = 0.189, p = 0.05). 

Discussion 

Species are present at a certain point because it fulfils three 
crucial conditions: (1) abiotic and (2) biotic conditions are 
favourable, and (3) the point is reachable for the species 
under consideration (Soberon & Peterson 2005). Condi¬ 
tion (2) refers to interactions with other species that mod¬ 
ify the species’ ability to maintain populations. This em¬ 
braces a wide range of co-occurring taxa that may serve 




Figure 2. Overlap between predicted distributions of the Pelomedusa species as determined by their ecological niche models. Scale 
bars = 1,000 km. 
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Figure 3. Major disjunct overpredicted areas for Pelomedusa 
barbata (black), P. kobe (grey), P. olivacea (yellow), and P. vari- 
abilis (orange). Red dots: records of Pelomedusa falling outside 
the modelled species’ ranges and outside the overpredicted ar¬ 
eas. Red lines enclose major overpredicted regions containing 
Pelomedusa records. For explanation of numbered areas see text. 
Scale bar = 1,000 km. 


as food, competitors, enemies or as vectors or triggers of 
diseases. The performed analyses in this study focused on 
abiotic conditions and information on the biotic interac¬ 
tion of each species is lacking. However, due to the fact 
that biotic factors robustly correlate with abiotic factors 
(Owens et al. 2012) and their impact disappears in large- 
scale analyses (Soberon & Nakamura 2009), we consider 
our developed models as meaningful approximations. 

The apparent correlation between the environments 
modelled as suitable for each species and biogeographical 
regions of sub-Saharan Africa indicates their congruency 
with general patterns. This suggests that the revealed eco¬ 
logical niches of the different taxa reflect real differences. 
Our statistical analyses (ANOVA, MANOVA, PDA) also 
corroborated the ecological niche models, supporting the 
notion that these species inhabit significantly different cli¬ 
matic envelopes (Fig. 4). Moreover, the fact that the relative 
contributions of each environmental variable and altitude 
are species-specific (Table r) suggests that the species are 
not ecologically compatible. Such ecological differentia¬ 
tion, together with allopatric or parapatric distribution of 
the species (Fig. r), indicates habitat isolation that prevents 
or limits gene flow. In spite of this, large areas of potential 
sympatry were revealed for a few species pairs (Fig. 2, Ta¬ 
ble 2). In line with our modelling results, the study by Pet- 
ZOLD et al. (2or4) already suggested the sympatric occur¬ 
rence of up to three species of Pelomedusa in the Arusha 
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Figure 4. Discriminant analysis of the climatic spaces of eight species of Pelomedusa. Percentages of the explained variation are given 
for both axes (CV - coefficient of variation). Total variation explained: 80.8%. Group centroids in red. 
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region of Tanzania, which was confirmed by our models 
(Fig. 2). We suggest that intensive surveys should be con¬ 
ducted in the potential regions of sympatry to confirm or 
reject sympatric or even syntopic occurrences of different 
Pelomedusa species there. 

We also identified disjunct overpredicted areas of consid¬ 
erable size in which unspecific Pelomedusa records are lo¬ 
cated (Fig. 3). Such a situation can be indicative of previously 
unrecognised species (Raxworthy et al. 2003, 2007), and 
the genetic and morphological screening of helmeted ter¬ 
rapins from these areas would be the consequent next step. 

However, there are also caveats regarding the reliability 
of our ecological niche models. Introduced populations of 
P. subrufa sensu stricto are widely distributed in northern, 
western, and southern Madagascar (cf the EMYSystem 
map and Petzold et al. 2014). Yet, our model for P. sub¬ 
rufa did not predict the occurrence of the species there 
(Fig. iB), but for other areas of Madagascar. This under¬ 
lines that ecological niche models are approximations that 
should not be trusted blindly. 

Even if our models are used under this premise for as¬ 
sessing the conservation status, they clearly indicate that 
the lUCN Threat Categories of the studied Pelomedusa 
species need to be revised. Until now, all helmeted terra¬ 
pins are still treated as a single widely distributed species 
by the lUCN, which is why helmeted turtles are not listed 
in the lUCN Red List of Threatened Species (lUCN 2014a). 
This corresponds to the Threat Category “Least Concern, 
LR/lc” (Boycott & Bourquin 2008). 

The extent of occurrence (EOO) and the area of occu¬ 
pancy (AOO) are both used as measures for assessing the 
conservation status of a species (Mace et al. 2008, lUCN 
2014b). According to our calculations, the EOOs and 
AOOs (Table 3) per se do not qualify any studied species 
for red-listing, because then EOOs < 20,000 km^ or AOOs 
< 2,000 km^ would be required. However, we assume that 
our calculations for the AOOs using the approach of Bach¬ 
man et al. (2011) are by far too optimistic. All Pelomedusa 
species depend on the presence of freshwater habitats and 
are therefore neither common nor widespread all over the 
inferred areas of occurrence. Moreover, surface freshwater 
is scarce and dwindling in the ranges of some Pelomedusa 
species. This is, for instance, true for P. barbata, a species 
endemic to the southwestern Arabian Peninsula. Petzold 
et al. (2014) already pointed out that this species is most 
probably endangered, and our AOO value does not neces¬ 
sarily contradict this view. 

As is obvious, our results constitute a first effort to¬ 
wards a better understanding of the distribution ranges of 
Pelomedusa species. The resulting models and statistical 
analyses are “under development” and awaiting new field 
data that could improve details and resolution. 
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Note added in proof 

In the meantime, Pelomedusa neumanni has been recorded from 
southern Ethiopia (Omo Region), i.e. from the immediate prox¬ 
imity of its predicted ecological niche. However, further to the 
north in Ethiopia, two deeply divergent genetic lineages of P. so- 
malica (two distinct candidate species) were found in regions 
within the ecological niche of P. neumanni (Oromia Region). At 
the eastern border of the predicted ecological niche of P. subrufa 
sensu stricto in South Africa, this species has been found at an¬ 
other site (Mpumalanga Province), so that the latter species has 
now been recorded twice in northeastern South Africa (Fritz et 
al. 2015 ). 
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Supplementary figure Si. Minimum convex polygons enclosing 
records used to develop the ecological niche models. 
Supplementary figure S 2 . Records of Pelomedusa from ElerpNET 
(http://www.herpnet.org), EMYSystem (http://emys.geo.orst. 
edu/), and Petzold et al. ( 2014 ). 

Supplementary figures S 3 -S 10 . Environments modelled as suit¬ 
able for Pelomedusa olivacea, P. neumanni, P. galeata, P. barbata, 
P. schweinfurthi, P. subrufa, P. variabilis, and P. kobe using MAX- 
ENT. 

Supplementary table Si. Additional records of Pelomedusa species 
used for ENM. 

Supplementary table S 2 . Statistics of univariate one-way analysis 
of variance (ANOVA) among Pelomedusa species for each envi¬ 
ronmental variable and altitude. 
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• Pelomedusa barbate 

• Pelomedusa galeata 

• Pelomedusa kobe 

• Pelomedusa neumanni 

• Pelomedusa ollvacea 

• Pelomedusa schweinfurthi 

• Pelomedusa subrufa 

• Pelomedusa variabilis 


Supplementary figure SI. Minimum convex polygons enclosing records used to develop the ecological niche models. 
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Supplementary figure S2. Records for Pelomedusa from HerpNET (http://www.herpnet.org), EMYSystem (http://emys.geo.orst.edu/), 
and Petzold et al. (2014). The allochthonous populations of Pelomedusa subrufa sensu stricto on Madagascar were excluded. The 
internet resources were accessed on 16 April 2013. 
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Supplementary figure S3. Environments modelled as suitable for Pelomedusa olivacea using MAXENT. Minimum training presence 
= 0.074. 
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Supplementary figure S4. Environments modelled as suitable for Pelomedusa neumanni using MAXENT. Minimum training pres¬ 
ence = 0.162. 
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Supplementary figure S5. Environments modelled as suitable for Pelomedusa galeata using MAXENT. Minimum training presence 
= 0.017. 
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Supplementary figure S6. Environments modelled as suitable for Pelomedusa barbata using MAXENT. Minimum training presence 
= 0.351. 
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Supplementary figure S7. Environments modelled as suitable for Pelomedusa schweinfurthi using MAXENT. Minimum training 
presence = 0.297. 
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Supplementary figure S8. Environments modelled as suitable for Pelomedusa subrufa using MAXENT. Minimum training presence 
= 0.381. 
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Online Supplementary data - Distribution modelling and conservation assessment for Pelomedusa spp. 



Supplementary figure S9. Environments modelled as suitable for Pelomedusa variabilis using MAXENT. Minimum training presence 
= 0.152. 
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Online Supplementary data - Mario Vargas-RamIrez et al. 



Supplementary figure SIO. Environments modelled as suitable for Pelomedusa kobe using MAXENT. Minimum training presence 
= 0 . 280 . 
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Online Supplementary data - Distribution modelling and conservation assessment for Pelomedusa spp. 
Supplementary table SI. Additional records of Pelomedusa species used for ENM (data from HerpNET and EMYSystem). 


Species 

Provenance 

Source 

Pelomedusa barbata 

Saudi Arabia: N17°01.00 E42°50.00 

HerpNet 

Pelomedusa barbata 

Saudi Arabia: N17°01.20, E42°49.80 

EMYSystem 

Pelomedusa barbata 

Saudi Arabia: N18°00.00, E42°34.80 

EMYSystem 

Pelomedusa barbata 

Yemen: N13°01.20, E44°54.00 

EMYSystem 

Pelomedusa barbata 

Yemen: N13°06.00, E045°22.80 

EMYSystem 

Pelomedusa barbata 

Yemen: N13°34.80, E44°01.80 

EMYSystem 

Pelomedusa barbata 

Yemen: N15°30.00, E43°25.20 

EMYSystem 

Pelomedusa kobe 

Tanzania: S06°09.00, E35°34.20 

EMYSystem 

Pelomedusa kobe 

Tanzania: S4°52.80, E34°09.00 

EMYSystem 

Pelomedusa kobe 

Tanzania: S5°52.20, E34°55.20 

EMYSystem 

Pelomedusa kobe 

Tanzania: S6°10.20, E35°04.20 

HerpNet 

Pelomedusa kobe 

Tanzania: S6°10.20, E35°40.20 

EMYSystem 

Pelomedusa kobe 

Tanzania: S6°10.99, E35°45.00 

EMYSystem 

Pelomedusa kobe 

Tanzania: S6°27.00, E35°34.20 

EMYSystem 

Pelomedusa neumanni 

Kenya: N0°03.00, E34°46.80 

EMYSystem 

Pelomedusa neumanni 

Kenya: N0°12.00, E35°28.01 

HerpNet 

Pelomedusa neumanni 

Kenya: N0°12.00, E35°28.20 

EMYSystem 

Pelomedusa neumanni 

Kenya: N0°40.80, S35°07.20 

EMYSystem 

Pelomedusa neumanni 

Kenya: N0°42.00, E35°22.80 

EMYSystem 

Pelomedusa neumanni 

Kenya: Sl°09.00, E35°46.20 

EMYSystem 

Pelomedusa neumanni 

Tanzania: S2°20.00, E34°50.00 

HerpNet 

Pelomedusa neumanni 

Tanzania: S2°21.31, E35°12.57 

HerpNet 

Pelomedusa neumanni 

Tanzania: S3°21.24, E35°36.73 

EMYSystem 

Pelomedusa olivacea 

Burkina Faso: N12°57.00; Wl°04.98 

EMYSystem 

Pelomedusa olivacea 

Burkina Faso: N13°51.00, W2°12.00 

EMYSystem 

Pelomedusa olivacea 

Nigeria: N12°12.00, E6°40.20 

EMYSystem 

Pelomedusa olivacea 

Nigeria: N12°51.28, E7°42.68 

HerpNet 

Pelomedusa olivacea 

Nigeria: N13°48.00, E5°22.80 

EMYSystem 

Pelomedusa olivacea 

Nigeria: N14°09.01, E8°30.99 

HerpNet 

Pelomedusa olivacea 

Togo: N10°54.00, E0°24.00 

EMYSystem 

Pelomedusa schweinfurthi 

Democratic Republic of the Congo: N3°37.20, E28°34.20 

EMYSystem 

Pelomedusa schweinfurthi 

Democratic Republic of the Congo: N3°52.80, E29°12.00 

EMYSystem 

Pelomedusa schweinfurthi 

Democratic Republic of the Congo: N4°10.16, E29°30.41 

EMYSystem 

Pelomedusa schweinfurthi 

Democratic Republic of the Congo: N4°22.20, E29°18.00 

EMYSystem 

Pelomedusa schweinfurthi 

South Sudan: N4°24.00, E32°34.20 

EMYSystem 

Pelomedusa schweinfurthi 

South Sudan: N4°24.29, E32°32.83 

HerpNet 

Pelomedusa schweinfurthi 

South Sudan: N4°42.00, E31°55.20 

EMYSystem 

Pelomedusa variabilis 

Ghana: N5°39.00, E0°10.80 

EMYSystem 

Pelomedusa variabilis 

Ghana: N5°39.00, E0°11.00 

HerpNet 

Pelomedusa variabilis 

Ghana: N6°06.00, E0°01.02 

HerpNet 

Pelomedusa variabilis 

Ghana: N6°06.00, E0°01.20 

EMYSystem 

Pelomedusa variabilis 

Ghana: N8°46.77, Wl°28.082 

HerpNet 
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Online Supplementary data - Mario Vargas-RamIrez et al. 


Supplementary table S2. Statistics of univariate one-way analysis of variance (ANOVA) among Pelomedusa species for each environ¬ 
mental variable and altitude. AMT - annual mean temperature; MDR - mean diurnal range; IT - isothermality; MTWM - maximum 
temperature of warmest month; MTCM - minimum temperature of coldest month; AP - annual precipitation; PWQ - precipitation 
of wettest quarter; PDQ - precipitation of driest quarter; PWaQ - precipitation of warmest quarter; PCQ - precipitation of coldest 
quarter; ALT - altitude. 


Bioclimatic variable 

Mean Square 

F 

df 

Sig. 

AMT 

25810.27 

26.114 

7 

0.0001 

MDR 

4487.94 

9.618 

7 

0.0001 

IT 

1367.46 

102.64 

7 

0.0001 

MTWM 

22482.36 

16.082 

7 

0.0001 

MTCM 

57239.33 

52.854 

7 

0.0001 

AP 

1553605.7 

14.792 

7 

0.0001 

PWQ 

291588.09 

16.017 

7 

0.0001 

PDQ 

25456.49 

12.253 

7 

0.0001 

PWaQ 

100033.53 

7.904 

7 

0.0001 

PCQ 

349107.62 

20.107 

7 

0.0001 

ALT 

4274818.72 

11.782 

7 

0.0001 
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